__device__ int RayBoxIntersect(float* tMin, float* tMax, float eye[3], float dir[3], float box0[3], float boxSize[3], int firstSeg[3], int lastSeg[3])
{
	float tmin=-FLT_MAX;
	float tmax=FLT_MAX;
	for(int i=0; i<3; i++)
	{
		if(dir[i]>ZERO_EPS)
		{
			tmin = max(tmin, (box0[i]+firstSeg[i]*boxSize[i]-eye[i])/dir[i]);
			tmax = min(tmax, (box0[i]+ lastSeg[i]*boxSize[i]-eye[i])/dir[i]);
		}
		else if(dir[i]<-ZERO_EPS)
		{
			tmin = max(tmin, (box0[i]+ lastSeg[i]*boxSize[i]-eye[i])/dir[i]);
			tmax = min(tmax, (box0[i]+firstSeg[i]*boxSize[i]-eye[i])/dir[i]);
		}
		if(tmin>tmax||tmax<0)	return 0;
	}

	(*tMin) = tmin;
	(*tMax) = tmax;
	return 1;

}

	//Normalize the coeffs
__device__ void normalizeCoeff(float* c, int n)
{
	float maxC=fabs(c[0]);
	for(int i=1; i<=n; i++)
	{
		maxC = max(maxC, fabs(c[i]));
	}
	maxC = 1.0/maxC;
	for(int i=0; i<=n; i++)
		c[i] *= maxC;
}

#if _KI_BEZ_
__device__ void deCasteljauSub(float u, int n, float* ori, float* L, float* R)
{
	float u1 = 1-u;
	int o=0, k1;
	float nn[MAX_DEG_INTERPOLATION+1][MAX_DEG_INTERPOLATION+1];


	for (int k=0; k<=n; k++)
		nn[0][k] = ori[o++];
	for (int k=1; k<=n; k++)
	{
		k1=k-1;
		for(int l=0; l<=n-k; l++)
			nn[k][l] = u1*nn[k1][l] + u*nn[k1][l+1];
	}
	for (int k=0; k<=n; k++)
	{
		L[k] = nn[k][0];
		R[k] = nn[n-k][k];
	}
}



__device__ int KnotInsertBezAll(float* c, int n, float* root)
{
	//ref: Real-Time Rendering of Trimmed Surfaces
	//1989 ACM SIGGRAPH by Alyn Rockwood
	//see page 110 one the right column "Root Finder Algorithm"

	//Normalize the coeffs
	normalizeCoeff(c, n);

	int rootN =0;
	//if c(0)=0 return t(1);
	if (fabsf(c[0])<cRenderParam.funcErr)
	{
		(root[rootN++]) = 0;
	}
	float lastR = 0.0f;
	if(fabsf(c[n])<cRenderParam.funcErr)
		lastR = 1.0f;


	int k=0;
	float L[MAX_DEG_INTERPOLATION+1];
	float R[MAX_DEG_INTERPOLATION+1];
	float actR, curR, scale, oldCR, eps=1.0e-4, delta;
	bool onL;
	float *LR;

	for (int i=0; i<=n; i++)
	{
		L[i] = c[i];
		R[i] = 1;
	}
	actR = curR = scale = oldCR = 1;
	while(true)
	{
		if (fabs(R[0])<cRenderParam.funcErr)
		{
			(root[rootN++]) = actR;

			actR += cRenderParam.deltaEps;
			deCasteljauSub(actR+cRenderParam.deltaEps, n, c, L, R);
			curR = oldCR = actR;
			scale = 1;

			if(R[0]*R[1]<0)
			{
				R[0] = -R[0];
			}
			for (int i=0; i<=n; i++)
			{
				L[i] = 1;
			}
		}
		
		for(k=0; k<n; k++)
		{
			if (L[k]*L[k+1]<0)	break;
		}
		if (k<n)
		{
			//check the left side.
			onL = true;
			LR = L;
		}
		else
		{
			for(k=0; k<n; k++)
			{
				if (R[k]*R[k+1]<0)	break;
			}
			if (k<n)
			{
				//check the right side.
				onL = false;
				LR = R;
			}
			else
			{
				if(lastR==1.0f)
				{
					(root[rootN++]) = 1;
				}
				return rootN;
			}
		}
		delta = LR[k]/((LR[k]-LR[k+1])*n);
		delta = max(delta, eps);
		curR = (k+0.0)/n + delta;
		deCasteljauSub(curR, n, LR, L, R);
		if (onL)
		{
			scale *= oldCR;
			actR -= scale*(1-curR);
		}
		else
		{
			scale *= (1-oldCR);
			actR += scale*curR;
		}
		oldCR = curR;
	}
}


__device__ float KnotInsertBez(float* c, int n)
{
	//ref: Real-Time Rendering of Trimmed Surfaces
	//1989 ACM SIGGRAPH by Alyn Rockwood
	//see page 110 one the right column "Root Finder Algorithm"

	//Normalize the coeffs
	normalizeCoeff(c, n);

	//if c(0)=0 return t(1);
	if (fabsf(c[0])<cRenderParam.funcErr)	return 0;


	int k=0;
	float L[MAX_DEG_INTERPOLATION+1];
	float R[MAX_DEG_INTERPOLATION+1];
	float actR = 0, eps=cRenderParam.deltaEps/*1.0e-4*/;
	float curR, scale, oldCR, delta;
	bool onL;
	float *LR;

	for (int i=0; i<=n; i++)
	{
		L[i] = c[i];
		R[i] = 1;
	}
	curR = scale = oldCR = 1;

	//int its = 0;
	while(fabs(R[0])>cRenderParam.funcErr/*&&its<cRenderParam.funcErr.nIter*/)
	{	
		for(k=0; k<n; k++)
		{
			if (L[k]*L[k+1]<0)	break;
		}
		if (k<n)
		{
			//check the left side.
			onL = true;
			LR = L;
		}
		else
		{
			for(k=0; k<n; k++)
			{
				if (R[k]*R[k+1]<0)	break;
			}
			if (k<n)
			{
				//check the right side.
				onL = false;
				LR = R;
			}
			else
			{
				return -1;
			}
		}
		delta = LR[k]/((LR[k]-LR[k+1])*n);
		delta = max(delta, eps);
		curR = (k+0.0)/n + delta;
		if (onL)
		{
			curR *= oldCR;
			deCasteljauSub(curR, n, c, L, R);
			oldCR = curR;
		}
		else
		{
			actR += scale*oldCR;
			scale *= (1-oldCR);
			for(int i=0; i<=n; i++)
			{
				c[i] = R[i];
			}
			deCasteljauSub(curR, n, R, L, R);
			oldCR = curR;
		}
		//its++;
	}
	actR += scale*oldCR;
	return actR;
}

#else
__device__ float knotAverage(float* t, int deg, int k)
{
	float aver = 0.0f;
	for (int i=k+1; i<=k+deg; i++)		aver += t[i];
	return	aver/deg;
}
//__device__ bool converged(float* t, int deg, int k, float* xList, int it)
//{
//	float eps = xList[it]-xList[it-1];
//	//eps /= (t[k+deg]-t[k]);
//	return (fabs(eps)<cRenderParam.funcErr);
//}

__device__ int KnotInsertAll(float* c, int d, float* root)
{
	//Normalize the coeffs
	normalizeCoeff(c, d);

 	float t[MAX_KNOT_ALL];

	int rootN = 0;
	//if c(0)=0 return t(1);
	if (fabsf(c[0])<cRenderParam.funcErr)
	{
		(root[rootN++]) = 0;
	}
	float lastR = 0.0f;
	if(fabsf(c[d])<cRenderParam.funcErr)
		lastR = 1.0f;


	//Bezier knot sequence .
	int n = 2*(d+1);
	for (int i=0;  i<=d; i++)		t[i] = 0;
	for (int i=d+1; i<n; i++)		t[i] = 1;

	n = d+1;
	int k=1, mu, it=0;
	float xPre, x, alpha;
	for (; it<MAX_KNOT_ITER_ALL; it++)
	{
		while((k<n)&&(c[k-1]*c[k]>0))	k++;
		if (k==n)			break;

		x = knotAverage(t, d, k) - c[k] * (t[k+d]-t[k])/((c[k]-c[k-1])*d);
		//xList[it] = x;

		if (it>=d&&fabsf(x-xPre)<cRenderParam.deltaEps)
		{
			(root[rootN++]) = x;
			k++;
			continue;
		}

		//Identify the knot interval [t_i t_(i+1)] where x lies in. 
		mu = k;
		while(x>=t[mu+1])	mu++;
		for (int i=n; i>=mu+1; i--)	c[i] = c[i-1];
		for (int i=mu; i>=mu-d+1; i--)
		{
			alpha = (x-t[i])/(t[i+d]-t[i]);
			c[i] = (1-alpha)*c[i-1] + alpha*c[i];
		}

		//Insert x to the knot interval;
		for (int i=n+d+1; i>mu+1; i--)	t[i] = t[i-1];
		t[mu+1] = x;
		xPre = x;
		n++;
	}

	if(MAX_KNOT_ITER_ALL==it&&k!=n)
	{
		(root[rootN++]) = x;
	}
	if(k!=n-1&&lastR==1.0f)
	{
		(root[rootN++]) = 1;
	}
	return rootN;
}

__device__  float KnotInsert(float* c, int d)
{
	//Normalize the coeffs
	normalizeCoeff(c, d);

	//if c(0)=0 return t(1);
	if (fabsf(c[0])<cRenderParam.funcErr)	return 0;

	//Bezier knot sequence .
 	float t[MAX_KNOT];

	int n = 2*(d+1);
	for (int i=0;  i<=d; i++)		t[i] = 0;
	for (int i=d+1; i<n; i++)		t[i] = 1;

	n = d+1;
	int k=1, mu;
	float xPre, x, alpha;
	for (int it=0; it<MAX_KNOT_ITER; it++)
	{
		while((k<n)&&(c[k-1]*c[k]>0))	k++;
		if (k==n)
		{
			return -1;
		}

		x = knotAverage(t, d, k) - c[k] * (t[k+d]-t[k])/((c[k]-c[k-1])*d);
		//xList[it] = x;

		if (it>=d&&fabsf(x-xPre)<cRenderParam.deltaEps)
		{
			return x;
		}

		//Identify the knot interval [t_i t_(i+1)] where x lies in. 
		mu = k;
		while(x>=t[mu+1])	mu++;
		for (int i=n; i>=mu+1; i--)	c[i] = c[i-1];
		for (int i=mu; i>=mu-d+1; i--)
		{
			alpha = (x-t[i])/(t[i+d]-t[i]);
			c[i] = (1-alpha)*c[i-1] + alpha*c[i];
		}

		//Insert x to the knot interval;
		for (int i=n+d+1; i>mu+1; i--)	t[i] = t[i-1];
		t[mu+1] = x;
		xPre = x;
		n++;
	}

	return x;
}
#endif



